%--------------------------------------------------------------------------
% Model Equations
%--------------------------------------------------------------------------

external_function(name = int_capital,           nargs = 10);
external_function(name = int_money,             nargs = 14);
external_function(name = int_profit,            nargs = 14);
external_function(name = int_liq,               nargs = 16);
external_function(name = int_acquisition,       nargs = 14);
external_function(name = int_partial,           nargs = 16);
external_function(name = int_surplus,           nargs = 16);
external_function(name = int_CM_reallocation,   nargs = 14);

external_function(name = int_2nd_price,         nargs = 12);
external_function(name = int_prod_disp,         nargs = 13);
external_function(name = average_log_normal,    nargs = 3);

model;

% --------------------------- Calculate the integrals and terms used in FOCs


q_capital = int_capital(L, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, CHI_PI, CHI_PI_k);

q_profit  = int_profit(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), CHI_PI, CHI_PI_k, tau_k);

q_money   = int_money(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), CHI_PI, CHI_PI_k, tau_k);

q_liq     = int_liq(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), CHI_k, CHI_PI, CHI_PI_k, tau_k, Z / K(-1));


% -------------------------- FOCs for capital, money, labor, consumption, and entry
1 / w   = BETA / w(+1) * (1 - tau_k) * B(+1) * (q_average_prod
 + ALPHA_0(+1) * ALPHA * (1 - THETA) * q_capital(+1)
 + ALPHA_0(+1) * ALPHA * THETA * CHI_k * (1 - DELTA) * q_money(+1)
 + ALPHA_0(+1) * ALPHA * THETA * CHI_PI_k * q_profit(+1)
 ) + BETA * (1 - DELTA) / w(+1);

%%%% monetary/non-monetary equilibriium
@#if case == 6
    %(1 + Mu) / w * Z  = BETA / w(+1) * Z(+1) * (1 + (1 - tau_k) * B(+1) * ALPHA * THETA * q_money);
    1 / w  = BETA / w(+1) / (1 + pi(+1)) * (1 + iota(+1));
@#else
    Z = 0;
@#endif

w       = (B / A / (1 - ETA))^((ETA - 1) / ETA) * ETA;

c_star^(-SIGMA) = XI / (1 - tau_h) / w;

surplus_net * w / XI = K(-1) * int_surplus(Z / K(-1), 1, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA,  DELTA, Chi_q(-1), CHI_PI, CHI_PI_k, CHI_k, tau_k) 
- iota * Z;

%%%% Entry condition (with log-Normal distribution of cost) 
ALPHA_0 = 0.5 * (1 + erf((log(surplus_net) - c_mu) / c_sig / sqrt(2)));

IM = average_log_normal(surplus_net, c_mu, c_sig);


% -------------------------- Monetary Policy

@#if case == 6 

    ALPHA_0 * Z  * (1 + pi) = ALPHA_0(-1) * Z(-1)  * (1 + Mu(-1));

    (1 + iota(-1)) = (1 + r) * (1 + pi);

   % (1 + iota) * BETA / (1 + pi(+1)) / w(+1) = 1 / w;
    
    iota = (1 - tau_k) * B * ALPHA * THETA * q_money;


@#else

    pi = Mu;

    r = 1 / BETA - 1; 

    iota = (1 + r) * (1 + pi) - 1; % could be any number, does not matter when money does not circulate

@#endif

% -------------------------- Market Clearing

L * (1 - tau_k) * B = Z / K(-1) + CHI_k * (1 - DELTA) - (1 - DELTA) * (1 - Chi_q(-1));

E_p = ALPHA_0 * ALPHA * K(-1) 
    * int_partial(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA,  DELTA, Chi_q(-1), CHI_k, CHI_PI, CHI_PI_k, tau_k, Z / K(-1));

E_a =  ALPHA_0 * ALPHA * K(-1) * int_acquisition(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), CHI_PI, CHI_PI_k, tau_k);
    
R_ratio * (E_a + E_p  + INV) = (E_a + E_p);

P_share * (E_a + E_p) = E_p ; 

K_bar = q_average_prod * K(-1) + ALPHA_0 *  ALPHA * K(-1) * q_capital + ALPHA_0 *  ALPHA * K(-1) * q_liq;                                   

Y   = c_star + G + K - (1 - DELTA) * K(-1);

Y   = B * K_bar / (1 - ETA);

INV = K - (1 - DELTA) * K(-1);

H   = (ETA / w)^(1 / (1 - ETA)) * A * K_bar;

Productivity * A^(1 - ETA) * K(-1)^(1 - ETA) * H^ETA = Y;

TFP  = Productivity * A^(1 - ETA);

DM_price    = int_2nd_price(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), tau_k);

prod_disp   = int_prod_disp(L, B, a, a, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, DELTA, Chi_q(-1), tau_k, ALPHA_0 * ALPHA);

% -------------------------- Observations
G_w     = log(w) - log(STEADY_STATE(w));

G_pi    = log(1 + pi) - log(1 + STEADY_STATE(pi));

G_iota  = iota - STEADY_STATE(iota); 

G_R_ratio   = log(R_ratio) - log(STEADY_STATE(R_ratio));

G_P_share   = log(P_share) - log(STEADY_STATE(P_share));

G_prod  = log(Productivity) - log(STEADY_STATE(Productivity));

G_TFP   = log(TFP) - log(STEADY_STATE(TFP));

@#if case == 6
   G_Z     = log(Z) - log(STEADY_STATE(Z));

@#else
   G_Z     = 0;
@#endif


G_output     = log(Y + IM * w * ZETA - G) - log(STEADY_STATE(Y + IM * w * ZETA - G)); % because output in the data is private output

G_INV   = log(INV) - log(STEADY_STATE(INV));

G_C     = log(c_star + IM * w * ZETA) - log(STEADY_STATE(c_star + IM * w * ZETA));

G_H     = log(H + IM) - log(STEADY_STATE(H + IM));

G_DM_price = log(DM_price) - log(STEADY_STATE(DM_price));

G_prod_disp  = log(prod_disp) - log(STEADY_STATE(prod_disp));

% -------------------------- Shocks

log(Chi_q) - log(CHI_q) = RHO_CREDIT  * (log(Chi_q(-1)) - log(CHI_q)) + SIGMA_CREDIT * e_CREDIT;

log(1 + Mu)   - log(1 + MU)   = RHO_MU  * (log(1 + Mu(-1))  - log(1 + MU)) - SIGMA_MU * e_Mu;

log(A) - log(A_SS) = RHO_A * (log(A(-1)) - log(A_SS)) + SIGMA_A * e_A;

log(G) - log(G_SS) = RHO_G * (log(G(-1)) - log(G_SS)) + SIGMA_G * e_G;

end;
